import numpy as np


def create(n):
    """
    Create a n-dim Extended Powell singular function

    :param n: dimensionality
    :return: x0, f, g, G
    """
    if n % 4 != 0:
        raise ValueError("n must be a multiple of 4")

    x0 = np.zeros(n, dtype=np.float64)
    x0[0::4] = 3
    x0[1::4] = -1
    x0[2::4] = 0
    x0[3::4] = 1

    def f(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        x1 = x[0::4]
        x2 = x[1::4]
        x3 = x[2::4]
        x4 = x[3::4]
        f1 = (x1 + 10 * x2) ** 2
        f2 = 5 * (x3 - x4) ** 2
        f3 = (x2 - 2 * x3) ** 4
        f4 = 10 * (x1 - x4) ** 4
        return np.sum(f1 + f2 + f3 + f4)

    def g(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        grad = np.zeros(n, dtype=np.float64)
        x1, x2, x3, x4 = x[0::4], x[1::4], x[2::4], x[3::4]
        grad[0::4] = 2 * (x1 + 10 * x2) + 40 * (x1 - x4) ** 3
        grad[1::4] = 20 * (x1 + 10 * x2) + 4 * (x2 - 2 * x3) ** 3
        grad[2::4] = 10 * (x3 - x4) - 8 * (x2 - 2 * x3) ** 3
        grad[3::4] = -10 * (x3 - x4) - 40 * (x1 - x4) ** 3
        return grad

    def h(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        he = np.zeros((n, n), dtype=np.float64)
        o = 0
        while o < n:
            x1, x2, x3, x4 = x[o], x[o + 1], x[o + 2], x[o + 3]
            he[o + 0, o + 0] = 2 + 120 * (x1 - x4) ** 2
            he[o + 0, o + 1] = 20
            he[o + 0, o + 3] = -120 * (x1 - x4) ** 2
            he[o + 1, o + 0] = 20
            he[o + 1, o + 1] = 200 + 12 * (x2 - 2 * x3) ** 2
            he[o + 1, o + 2] = -24 * (x2 - 2 * x3) ** 2
            he[o + 2, o + 1] = -24 * (x2 - 2 * x3) ** 2
            he[o + 2, o + 2] = 10 + 48 * (x2 - 2 * x3) ** 2
            he[o + 2, o + 3] = -10
            he[o + 3, o + 0] = -120 * (x1 - x4) ** 2
            he[o + 3, o + 2] = -10
            he[o + 3, o + 3] = 10 + 120 * (x1 - x4) ** 2
            o += 4
        return he

    return x0, f, g, h
